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Abstract 

We present first-principles calculations of elastic properties of multilayered two-dimensional crys¬ 
tals such as graphene, h-BN and 2 J/-M 0 S 2 which shows that their Poisson’s ratios along out-of- 
plane direction are negative, near zero and positive, respectively, spanning all possibilities for sign 
of the ratios. While the in-plane Poisson’s ratios are all positive regardless of their disparate elec¬ 
tronic and structural properties, the characteristic interlayer interactions as well as layer stacking 
structures are shown to determine the sign of their out-of-plane ratios. Thorough investigation of 
elastic properties as a function of the number of layers for each system is also provided, highlighting 
their intertwined nature between elastic and electronic properties. 

PACS numbers: 62.20.de, 62.20.dj, 73.21.Ac 
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Under uniaxial stress, Poisson’s ratio defined by the ratio of the strain in the transverse 
direction (et) to that of the longitudinal direction (e^), ly = —et/ei, measures the fundamental 
mechanical responses of solids against external loads [l-5|. It has strong correlation with 
atomic packing density, atomic connectivity and structural phase transition M . The 
theory of elasticity allows values of Poisson’s ratio of an isotropic material ranging from —1 to 
0.5, i.e., from extremely compressible to incompressible materials [l, 5|. Thus, when a solid 
is subjected to a uniaxial compression, it expands (z/ > 0), remains to be the same [u = 0), 
and shrinks {u < 0) in the transverse direction depending on the sign of Poisson’s ratio. 
Typically, different Poisson’s ratio or its sign indicates dramatic variations in mechanical 
properties. For example, when isothermal modulus is extremely larger than shear modulus, 
the material reaches its incompressible limit as shown in most liquids or rubber [u ~ 0.5) 
and in the opposite case, re-entrant foams and related structures show the negative u or 
auxetic property js-O]. The Poisson’s ratio of common solid state crystals usually falls in 
the range of 0 < < 0.5 while gases and cork have i/ ~ 0 

Anisotropic materials with directional elastic properties often shows more dramatic vari- 


the experimental realization of graphene 
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ations in their Poisson’s ratios such as the directional auxetic property . In this regard. 


the thinnest and the strongest mate- 


15| . now offers a new platform to understand electronic and elastic properties of 


well-dehned anisotropic materials and their heterostructures. Even though the Young’s 
modulus and Poisson’s ratio of graphene have been studied quite thoroughly 12M2l|. those 
along the out-of-plane direction for its few-layered forms have barely been known. Neither 
do for all the other available two-dimensional crystals. Since electronic properties of layered 



calculated fundamental mechanical properties of three representative van der Waals (vdW) 
crystals along all crystallographic directions of their few-layered structures. 

In this work, we present a theoretical study using a hrst-principles approach on the 
elastic properties of layered two-dimensional crystals, including graphene, h-BN and 2H- 
M 0 S 2 , in which the vdW energy is one of the governing interactions between their layers 
while they exhibit very different electronic properties. We hnd that the Poisson’s ratios 


2 




































of graphene, h-BN and 2iir-MoS2 along out-of-plane direction are negative, near zero and 
positive, respectively, whereas their in-plane Poisson’s ratios are all positive. The diverseness 
of out-of-plane Poisson’s ratio is attributed to their disparate electronic properties as well 
as stacking structures. Thorough investigation on their elastic properties while varying the 
number of layers are also reported. 

We hrst consider graphene with AB stacking, h-BN and 2i7-MoS2 with AA' stacking. 
All the three have symmetry. Generally, for a material with symmetry, the stiffness 
tensor without shear part can be written with four independent parameters. 
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with the choice of 2 ; as the axis for the three-fold rotational symmetry [1[. Here, cxj and 
are the stress and strain respectively along the i-th axis. The components of stiffness tensor 
can be obtained by differentiating the total energy in terms of strain; A = d‘^E<^oi/de\ = 
d'^Etotldel, B = d'^E^otlde^dey, C = and D = d‘^Etotjdel. By taking the 

inverse of the stiffness tensor, one can get the compliance tensor. 
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The subscripts i and o represent in-plane and out-of-plane respectively. Ei and Eo are the 
Young’s moduli along the x{y) and 2 ; axis respectively. There are two out-of-plane Poisson’s 
ratios; Uq is the Poisson’s ratio along the z axis when the stress is applied along the x or 
y directions while Uq = i^oEq/E^ is the Poisson’s ratio along the x or y direction when the 
stress is applied along the 2 ; direction. Ui is the in-plane Poisson’s ratio along the x{y) axis 
when the stress is applied along the y{x) axis. 

Using a first-principles approach based on density-functional theory with plane wave basis 
set [2^, we calculate total energies, Etot{^x,^y,^z), of all systems at 5 x 5 x 5 grid points 


in the strain space of {ex,ey,ez). To obtain the accurate binding energy and interlayer 


distance including the vdW energy, we have used the revised version 
correlation functional method developed by Vydrov and van Voorhis 


27( 1 of the nonlocal 


for reproducing both values following results from more accurate methods 


that is successful 
2911. In order to 
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FIG. 1: Lattice structures of (a) AS-stacked graphene (b) A^'-gtacked /i-BN and (c) 2 LI-MoS 2 . 
The parameter d is the interlayer distance of each structure. In 2 iL-MoS 2 , d is the vertical distance 
between Mo atoms in adjacent layers and di is the vertical intralayer sulfur-to-sulfur distance. 

reduce spurious interaction between neighbouring supercells, a large vacuum over 68 A is 
introduced and relatively high energy cut-off above 100 Rydberg as well as dense k-point 
grids up to 29 x 29 are used to converge the results. 

For the total energy calculations with tensile strain on all systems in which the layers are 
stacked along the 2 ; axis, a primitive cell with unit vectors, Ui = (a, 6 , 0 ), U 2 = (—a, 6 , 0 ), 
U 3 = (0, 0, c) are used [Fig. [I]. Strain along x and y axis is defined by e^, = (a — ao)/ao and 
~ 6 o)/ 6 o, where oq and 60 are the lattice parameters of the equilibrium structure. 
Strain across the layers along the 2 ; axis is defined by = (d — do)/do where d and do are 
the interlayer distance and that of the equilibrium structure, respectively. The calculated 
lattice parameters of oq and do for infinitely stacked bulk systems are 2.47 A, 2.52 A, 3.22 



do depending on the number of layers are reflected in our calculations. A, B, C and D in 
Eq. (1) are calculated by interpolating Etot on the strain space and Young’s modulus and 
Poisson’s ratio from compliance tensor in Eq. (2). 

Figure [2] shows Young’s moduli and Poisson’s ratios for the three materials with vari¬ 
ous number of layers. In-plane elastic constants, Ei and z/*, are barely dependent on the 
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FIG. 2: Elasticity constants of graphene (empty rectangle), /i-BN (filled rectangle) and 2ff-MoS2 
(cross) as a function of number of layers upto ten layers, in-plane (a) Young’s modulus and (b) 
Poisson’s ratio and out-of-plane (c) Young’s modulus and (d) Poisson’s ratio. 


number of layers [Figs. |2][a) and (b)]. Under in-plane tensile stress, the hexagonal network 
of atoms is deformed for graphene and h-BN while, for 2hf-MoS2, sulfur-to-sulfur distance 
across the plane along the z axis within one layer can also be deformed. Furthermore, the 
hexagonal structure of graphene and h-BN with rigid a bond supported by vr bond is stiffer 
than that of 2iJ-MoS2. This makes 2hf-MoS2 more flexible to applied stress resulting in 
lower in-plane Young’s modulus [Fig. [2](a)] and larger in-plane Poisson’s ratio [Fig. [2][b)]. 
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Young’s moduli across layers, Eo, increase for all of the three materials with the increase 
of the number of layers reflecting the accumulation of long-range interlayer van der Waals 
interaction [Fig. [2]^c)]. 

Contrary to similar behaviours between in-plane elastic properties of the three materials, 
out-of-plane elastic properties between those differ qualitatively [Fig. 2(d)]. The most 
notable one is that multilayered graphene structures have out-of-plane Poisson’s ratio as 
negative as u = —0.09 [Fig. [2] (d)] with slight dependence on the layer number variations. 
Materials with axial negative Poisson’s ratio have been reported during the last few decades 
such as foams with re-entrant atomic structures js-O] and those with non-axial one are 
shown in some simple cubic metals 3J, [3^. The present case is for the axial negative ratio 
in a layered material where the vdW interaction governs the binding between layers without 
re-entrant structure. More interestingly, h-BN shows very small out-of-plane Poisson’s ratio 
near zero crossing from negative to positive values as number of layers increases whereas 
2Ff-MoS2 has positive Poisson’s ratio as shown in Fig. 2(d). So, the three layered crystals 
have qualitatively different Poisson’s ratios spanning all possibilities of their signs. 

To understand the qualitative difference in out-of-plane Poisson’s ratios of the three 
layered systems, we hrst decompose the binding energy of bilayer systems into repulsive 
and attractive parts. Figure [3](a) shows the binding energy curve between two layers of 
graphene, ^bind(d) = | where is the total energy of bilayer (single 

layer) graphene. Eund{d) is shown as solid (no strain) and open (8.1% equibiaxial nominal 
strain) circles as a function of the distance, d, between the layers. The calculated binding 
energy within d = 4 ~ 9 A is well described by i?bind(d) ~ d~^. The fitting curves of d~‘^, 
which reflect asymptotic vdW interaction, are drawn in as a solid (dashed) line without 
(with) strain. The difference between the total energy and vdW energy for each case is also 
plotted in the same pl^ two curves on top, representing purely repulsive characteristics 
called Pauli repulsion |36|, l37| . It does not £t to any single power of but is well fit by 
an exponentially decaying function, Eundid) ~ [exp(dVa2) - 1] \ with a = 1.37 A for the 
case without strain. We note that the d~^ dependence of vdW energy of bilayer graphene. 


which was recently reported 


is valid only at a distance larger than 9 A 39[ therefore 


not relevant near equilibrium distance considered here. 

Figure [3](a) indicates that both vdW attraction and Pauli repulsion on bilayer graphene 
are enhanced under tensile strain. However, noting that the Pauli repulsion energy showing 
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FIG. 3: Binding energy as a function of interlayer distance, d, for bilayer (a) graphene and (b) 
/i-BN. The htting curves for attractive and repulsive parts are drawn in for unstrained case (solid 
line) and ex = ey = 0.08 equibiaxial strain case (dashed line). The insets show magnified views 
near equilibrium. 


exponential increase with d is ninch stiffer than attractive vdW energy, equilibrinm interlayer 
distance is critically sensitive to the change of the former than the latter. Thus, the equilib¬ 
rium interlayer distance under the strain is mainly determined not by the strain-enhanced 
vdW attraction but by the enhancement of the repulsion. In graphene, electronic states 
pointing away from the layers are composed of linear combinations of pz orbitals of atoms 
called TT orbitals and form the vr band Since the occupied electrons of vr orbitals in 


adjacent layers expel each other from their overlap region 42|, the enhanced repulsion with 
external strain shown in Fig. [3] (a) may indicate the strain-induced spatial variation of vr 
electrons pushing the two layers away while the vdW interaction still keeps their binding. 

We hnd that the in-plane strain indeed elongates the spatial distribution of electron 
density away from the layer making the Pauli repulsion increase over the vdW attraction. 
We calculate the spatial distribution of density of vr band along the 2 ; axis in a single layer 
graphene; Pn{z) = dE f dxdypE{x, y, z). Here, p _ e ( x , y, z) is the local density of state for 
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the vr band only and is summed over the in-plane unit cell {x and y axis). Graphene is located 
at z = 0 and the Fermi energy of the neutral system is denoted by Sp. Our ab initio calcula¬ 
tion result for the maximum of pT^i^) decreases with tensile strain while its tail increases im¬ 
plying that TT orbital spreads out along the 2 ; axis with strain. Quantitatively, we calculate the 
density-weighted length of tt orbital along the z axis using = f dz \z\pt,{z)/ f dz Ptt{z) 
that gives = 0.673 A without strain. We find that the value of Lt, indeed increases by 
0.6 % as the equibiaxial strain increases by 2 %, thereby explaining the value of the negative 
Poisson’s ratio near —0.1 along out-of-plane direction. This elongation can be understood 
simply by considering overlaps between neighbourig atomic orbitals. For a charge neutral 
graphene, the spatial distribution of vr orbitals along the perpendicular direction to the layer 
is contracted compared to orbitals of an isolated carbon atom because of overlap between 
nearby pz orbitals. In-plain tensile strain returns the carbon atoms in graphene back to 
isolated one so that the tt orbitals should be elongated. 

A simple tight-binding (TB) picture can corroborate the elongation of spatial distribution 
of TT orbitals under strain. Consider the Bloch wave function within the TB approximation, 
(pA{B)ik,r) = — Ra{b))i where subscript A{B) represents the 

sublattice index, Pa{b){.^— Ra{b)) is the normalized Pz orbital of the carbon atom at Ra{b), 
while Ra(b) runs the positions of atoms in the A{B) sublattice [^. With the nearest 
neighbor hopping, t, and the overlap, s = {ipAi.'i^ — Ra)\^b{A — Ra — Sj)), the vr-orbital 


is given by A) = T 7 ^ 0 yi(^,A) -f- (j)B{k,r) where f{k) = ^^cl 5j points to 

the three nearest neighbors. Considering s -C 1, = A— j^^d^k ln{k), where ln{k) = 

~ ~ I/(^)I(A 2 ® “ i® density-weighted length of the vr orbital at k and 

Sbz is the area of the first Brillouin zone (BZ). Here, Ip^ = f d^r |z||y 9 ^(F)p is the length of 
an isolated pz orbital and h = J dz’r ip\{r — RA)\z\ipB{A— Ra — (5i). If the maximum overlap 
between nearest neighbor pz orbitals is at \z\ = zq and (p’^(f— Ra)(Pb(p — Ra — di) is trivial 
elsewhere, then 4 ~ so that C(A;) ~ Ip^ — \f{k)\{lp^ — Zo)s and that ~ Ip^ — 4 (Ai|) 
where Ih = /bz i® straightforward to hnd that Ih > 0 and dlh/d\6i\ < 0 . 

Therefore, the above simple formulation for L^r implies that the out-of-plane distance of tt- 
orbital is shorter than that of bare Pz orbital and that the applied tensile strain can increase 
its distance. 


Now, let us compare elastic properties of multilayered h-BN with graphene. A previous 
study |43| shows that the ionic interaction energy in the h-BN is negligible in determining 






the electrostatic repulsion as well as dispersion forces so that the interlayer distance is very 
similar to graphite regardless of apparent difference in the static polarizability between the 
two layered materials. Our analysis, however, shows the elastic properties can be quite 
different; amplitude of out-of-plane Poisson’s ratio of h-BN is nowhere close to that of 
multilayer graphene but order of magnitude smaller [Fig. [2][d)]. Figure |2](b) shows that 
under tensile strain, vdW interaction increases as in graphene while the Pauli repulsion 
barely changes. For a bilayer graphene with AB stacking, half of the carbon atoms of one 
layer are right on top of the carbon atoms of the other layer so that the tails of pz orbitals 
from two layers directly overlap with each other. For a bilayer h-BN with AA' stacking, 
however, fully-hlled Pz orbital of nitrogen atom from one layer is on top of the empty Pz 
orbital of boron from the other layer so that the interlayer Puali repulsion is not as sensitive 
to the slight change of the length of pz orbitals as in the case of graphene. 

The length of Pz orbital of a single-layer h-BN, in fact, does change due to the strain in 
the same manner as that of graphene. Therefore, a negative interlayer Poisson’s ratio should 
appear in h-BN as multilayer graphene if the interlayer alignment of pz orbitals follows that 
of graphene. This can be realized by changing the stacking structure of h-BN from AA' to 
AA. We have computed the elastic properties for this artificial bilayer structure and found 
that the out-of-plane Poisson ratio is —0.12, thus conhrming our theory. On the other hand, 
the out-of-plane Poisson’s ratio of bilayer 2iJ-MoS2 is mainly determined by the flattening 
of each layer under tensile strain that gives positive value. Our calculation shows that the 
change of interlayer distance, d, of 2id-MoS2 in Fig. [T][c) in response to a given in-plane 
tensile stress mainly comes from the change of di; Adi ~ (3/4)Ad. 

In conclusion, we have studied the elastic properties of multilayered two-dimensional 
crystals including graphene, h-BN, and 2dd-MoS2, with interlayer van der Waals interaction 
properly taken into account. In-plane elastic properties are found to be barely dependent 
on the number of layers for all three materials. Our analysis reveals that graphene is a very 
peculiar axial auxetic material when in-plain strain is applied. The mechanism is attributed 
to quantum mechanical origin rather than to structural one such as re-entrant foam. In 
contrast, the Poisson’s ratio of h-BN with AA' stacking is found to be nearly zero and that 
of M 0 S 2 is positive. 
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